Study design:
#clear all inputs
rm(list = ls())
#Check for necessary packages
#list.of.packages <- c("ggplot2",
# "tidyr",
# "dplyr",
# "Hmisc",
# "qqplotr",
# "ggthemes",
# "fabricatr",
# "gridExtra",
# "grid",
# "kableExtra",
# "sjPlot",
# "cowplot",
# "ggpubr",
# "patchwork",
# "magick", #so you can use the savekable function
# "webshot", #so you can use the savekable function
# "lme4",
# "RcppEigen")
#should install packages that you don't have
#new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
#if(length(new.packages)) install.packages(new.packages)
#Load the package that contains the Cox Propotional Hazard Function
library(ggplot2)
library(tidyr)
library(dplyr)
library(Hmisc)
library(qqplotr)
library(ggthemes)
library(fabricatr)
library(gridExtra)
library(grid)
library(kableExtra)
library(sjPlot)
library(cowplot)
library(ggpubr)
library(patchwork)
#This imports the data from a CSV file and creates the data frame PrTC
Adult_Final<-read.table("2021_Adult_Firefly_Survival_binned3.csv", header=TRUE, sep=",", dec=".",na.strings=".")
#This sets up the color palette using tableau20: https://jrnold.github.io/ggthemes/reference/tableau_color_pal.html
site_colors <- c("USA: Montour Co, Bucknell Natural Area" = "#4E79A7",
"USA: Union Co, Bucknell Farm" = "#F28E2B",
"USA: Union Co, Bucknell Ropes Course" = "#E15759")
Data Used:
Abbreviations used:
#generate distribution of raw data
z <- ggplot(Adult_Final, aes(x=ElytralLength)) +
geom_histogram() +
xlab("Elytron length (mm)") +
ylab("N fireflies") +
theme(plot.title = element_text(hjust = 0.5))
#generate log distribution plot
y <- ggplot(Adult_Final, aes(x=log(ElytralLength))) +
geom_histogram() +
xlab("log(Elytron length) (mm)") +
ylab("N fireflies") +
theme(plot.title = element_text(hjust = 0.5))
grid.arrange(z, y, nrow= 1)
#generate tibble
t<- Adult_Final %>%
group_by(Location) %>%
summarise(N=n(), mean_length=round(mean(ElytralLength),2), SD = round(sd(ElytralLength), 2))
#generate nice looking table
kbl(t, col.names = c("Location", "n", "Mean elytron length (mm)", "SD"), align = "lccc") %>%
kable_classic(full_width = F, html_font = "Cambria")
| Location | n | Mean elytron length (mm) | SD |
|---|---|---|---|
| USA: Montour Co, Bucknell Natural Area | 298 | 9.11 | 0.96 |
| USA: Union Co, Bucknell Farm | 70 | 8.90 | 1.08 |
| USA: Union Co, Bucknell Ropes Course | 470 | 8.94 | 0.95 |
#generate labels for sites that include sample size
site_labels = c(paste("BNA", " (n = ", t$N[1], ")", sep=""),
paste0("BUF", " (n = ", t$N[2], ")"),
paste0("FDBCC", " (n = ", t$N[3], ")"))
#generate text placeholder
a <- ggplot(Adult_Final, aes(x=ElytralLength, fill=Location)) +
geom_histogram(alpha= 0.5, position="identity") +
scale_fill_manual(name = "Collection Location", labels = site_labels, values = site_colors)
legend.a <- cowplot::get_legend(a)
#generate histogram
b <- ggplot(Adult_Final, aes(x=ElytralLength, fill=Location)) +
geom_histogram(alpha= 0.5, position="identity") +
xlab("Elytron length (mm)") +
ylab("N fireflies") +
theme(plot.title = element_text(hjust = 0.5), legend.position = "none") +
scale_fill_manual(name = "Location", labels = site_labels, values = site_colors)
#generate violin plot
c <- ggplot(Adult_Final, aes(x=Location, y=ElytralLength, fill=Location)) +
geom_violin(alpha = 0.9) +
theme_classic() +
ylab("Elytron length (mm)") +
theme(plot.title = element_text(hjust = 0.5), legend.position = "none", axis.text.x=element_blank(), axis.title.x=element_blank()) +
geom_boxplot(width=0.2, alpha= 0.5, fill = "white") +
scale_fill_manual(values = site_colors)
#generate boxplot
d <- ggplot(Adult_Final, aes(x=Location, y=ElytralLength, fill=Location)) +
geom_boxplot(alpha = 0.9) +
geom_jitter(size = 0.1)+
theme_classic() +
ylab("Elytron length (mm)") +
theme(plot.title = element_text(hjust = 0.5), legend.position = "none", axis.text.x=element_blank(), axis.title.x=element_blank()) +
scale_fill_manual(values = site_colors)
grid.arrange(legend.a,b,c,d, nrow=2)
#generate tibble
t <- Adult_Final %>%
group_by(Location, SeasonDay) %>%
summarise(N=n(), mean_length=round(mean(ElytralLength),2), SD = round(sd(ElytralLength),2)) %>%
arrange(SeasonDay)
#make table nice in kable
kbl(t, col.names = c("Location", "Season Day", "n", "Mean elytron length (mm)", "SD"), align = "lcccc") %>%
kable_classic(full_width = F, html_font = "Cambria")
| Location | Season Day | n | Mean elytron length (mm) | SD |
|---|---|---|---|---|
| USA: Montour Co, Bucknell Natural Area | 9 | 92 | 9.55 | 0.91 |
| USA: Union Co, Bucknell Farm | 13 | 47 | 8.76 | 1.09 |
| USA: Montour Co, Bucknell Natural Area | 15 | 120 | 8.71 | 0.98 |
| USA: Union Co, Bucknell Ropes Course | 17 | 231 | 8.96 | 0.94 |
| USA: Montour Co, Bucknell Natural Area | 22 | 86 | 9.21 | 0.73 |
| USA: Union Co, Bucknell Ropes Course | 23 | 95 | 8.93 | 1.01 |
| USA: Union Co, Bucknell Farm | 30 | 23 | 9.21 | 1.01 |
| USA: Union Co, Bucknell Ropes Course | 31 | 90 | 8.94 | 0.93 |
| USA: Union Co, Bucknell Ropes Course | 36 | 54 | 8.90 | 0.97 |
#generate labels for seasondays that include sample size
t$SeasonDayLabel = paste(t$SeasonDay, " (", t$N, ")", sep ="")
seasonday_label <- c(
'9' = t$SeasonDayLabel[1],
'13' = t$SeasonDayLabel[2],
'15' = t$SeasonDayLabel[3],
'17' = t$SeasonDayLabel[4],
'22' = t$SeasonDayLabel[5],
'23' = t$SeasonDayLabel[6],
'30' = t$SeasonDayLabel[7],
'31' = t$SeasonDayLabel[8],
'36' = t$SeasonDayLabel[9]
)
#generate legend
e <- ggplot(Adult_Final, aes(x=ElytralLength, fill=Location)) +
geom_histogram() +
xlab("Elytron length (mm)") +
ylab("N fireflies") +
theme(plot.title = element_text(hjust = 0.5)) +
facet_wrap(~SeasonDay, labeller = as_labeller(seasonday_label)) +
scale_fill_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759"))
legend.e <- cowplot::get_legend(e)
#generate histogram
f <- ggplot(Adult_Final, aes(x=log(ElytralLength), fill=Location)) +
geom_histogram() +
xlab("Elytron length (mm)") +
ylab("N fireflies") +
ggtitle("SeasonDay (n)") +
theme(plot.title = element_text(hjust = 0.5), legend.position = "none") +
facet_wrap(~SeasonDay, labeller = as_labeller(seasonday_label), ncol=5) +
scale_fill_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759"))
#generate violin plot by seasonday
g <- ggplot(Adult_Final, aes(x=as.factor(SeasonDay), y=ElytralLength, fill = Location)) +
geom_violin(alpha = 0.9) +
theme_classic() +
ylab("Elytron length (mm)") +
xlab("Season Day") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90), legend.position = "none") +
geom_boxplot(width=0.2, fill = "white") +
scale_fill_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759"))
#generate violin plot by site
h <- ggplot(Adult_Final, aes(x=as.factor(SeasonDay), y=ElytralLength, fill = Location)) +
geom_violin(alpha = 0.9) +
theme_classic() +
ylab("Elytron length (mm)") +
xlab("Season Day") +
theme(plot.title = element_text(hjust = 0.5), legend.position = "none", axis.text.x = element_text(angle = 90)) +
facet_wrap(~Location, labeller = labeller(Location =
c("USA: Montour Co, Bucknell Natural Area" = "BNA",
"USA: Union Co, Bucknell Farm" = "BUF",
"USA: Union Co, Bucknell Ropes Course" = "FDBCC"))) +
geom_boxplot(width=0.2, fill = "white") +
scale_fill_manual(values = c("#4E79A7", "#F28E2B", "#E15759"))
lay1 <- rbind(c(2,2,2,1))
lay2 <- rbind(c(1,1,2,2,2))
g1 <- arrangeGrob(legend.e, f, layout_matrix = lay1)
g2 <- arrangeGrob(g,h, layout_matrix = lay2)
grid.arrange(g1, g2, nrow=2)
Testing the raw data:
#qq plots to assess normality
q <- ggplot(Adult_Final, aes(sample=ElytralLength, color=Location)) +
stat_qq_band(bandType = "pointwise", fill = "#8DA0CB", alpha = 0.4) +
stat_qq_line(color = "#8DA0CB") +
stat_qq_point() +
ggtitle("Season Day") +
facet_wrap(~SeasonDay) +
scale_color_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759")) +
theme(legend.position = "none")
#statistically test
shapiro.table <- NULL
for (day in levels(factor(Adult_Final$SeasonDay))){
#subset table
table.sub <- Adult_Final[Adult_Final$SeasonDay == day,]
test.result <- (shapiro.test(table.sub$ElytralLength))
info <- data.frame(SeasonDay = day, W = round(test.result[1]$statistic, 2), p = round(test.result[2]$p.value, 2))
shapiro.table <- rbind(shapiro.table, info)
}
#create kable filename
kable_filename = paste("Elytron_length_shapiro", ".png", sep="")
#generate nice table with kable
kbl(shapiro.table, row.names = FALSE, col.names = c("Season Day", "W", "P"), align = "ccc") %>%
kable_classic(full_width = F, html_font = "Cambria") %>%
save_kable(file = kable_filename, bs_theme = "cerulean", zoom = 7)
#make big plot with table
#read in kable table as image
r <- ggdraw() +
draw_image(kable_filename)
#stack the legend on top of the table
g3 <- arrangeGrob(legend.e, r, nrow =2)
#make the column width layout
shapiro_layout <- rbind(c(1,1,1,2))
#plot the graph next to the legend and table
grid.arrange(q, g3, nrow = 1, layout_matrix = shapiro_layout)
Conclusions:
This suggests that log transformation is not necessary to satisfy statistical assumptions of a linear relationship between the response variable (elytron length) and explanatory variables (season day, location) when looking at how these variables are related. However, log transformation is the standard when comparing body size metrics due to allometric scaling expectations. Thus, we will proceed with testing whether log transformation of elytron length alters expected linear relationships among variables (specifically body mass and survival). The fits may be no different or even better.
Testing log-transformed data:
#qq plots to assess normality
q <- ggplot(Adult_Final, aes(sample=log(ElytralLength), color=Location)) +
stat_qq_band(bandType = "pointwise", fill = "#8DA0CB", alpha = 0.4) +
stat_qq_line(color = "#8DA0CB") +
stat_qq_point() +
ggtitle("Season Day") +
facet_wrap(~SeasonDay) +
scale_color_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759")) +
theme(legend.position = "none")
#statistically test
shapiro.table <- NULL
for (day in levels(factor(Adult_Final$SeasonDay))){
#subset table
table.sub <- Adult_Final[Adult_Final$SeasonDay == day,]
test.result <- (shapiro.test(log(table.sub$ElytralLength)))
info <- data.frame(SeasonDay = day, W = round(test.result[1]$statistic, 2), p = round(test.result[2]$p.value, 2))
shapiro.table <- rbind(shapiro.table, info)
}
#create kable filename
kable_filename = paste("log_Elytron_length_shapiro", ".png", sep="")
#generate nice table with kable
kbl(shapiro.table, row.names = FALSE, col.names = c("Season Day", "W", "P"), align = "ccc") %>%
kable_classic(full_width = F, html_font = "Cambria") %>%
save_kable(file = kable_filename, bs_theme = "cerulean", zoom = 7)
#make big plot with table
#read in kable table as image
r <- ggdraw() +
draw_image(kable_filename)
#stack the legend on top of the table
g3 <- arrangeGrob(legend.e, r, nrow =2)
#make the column width layout
shapiro_layout <- rbind(c(1,1,1,2))
#plot the graph next to the legend and table
grid.arrange(q, g3, nrow = 1, layout_matrix = shapiro_layout)
Conclusions:
Elytral<-lm(ElytralLength~SeasonDay*Location, data=Adult_Final)
summary(Elytral)
##
## Call:
## lm(formula = ElytralLength ~ SeasonDay * Location, data = Adult_Final)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.98747 -0.62992 0.00396 0.64271 2.91253
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 9.48258 0.17763
## SeasonDay -0.02430 0.01112
## LocationUSA: Union Co, Bucknell Farm -1.06859 0.34167
## LocationUSA: Union Co, Bucknell Ropes Course -0.48694 0.23421
## SeasonDay:LocationUSA: Union Co, Bucknell Farm 0.05070 0.01821
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course 0.02206 0.01279
## t value Pr(>|t|)
## (Intercept) 53.385 < 2e-16 ***
## SeasonDay -2.186 0.02910 *
## LocationUSA: Union Co, Bucknell Farm -3.128 0.00182 **
## LocationUSA: Union Co, Bucknell Ropes Course -2.079 0.03792 *
## SeasonDay:LocationUSA: Union Co, Bucknell Farm 2.784 0.00550 **
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course 1.724 0.08505 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9639 on 832 degrees of freedom
## Multiple R-squared: 0.01736, Adjusted R-squared: 0.01145
## F-statistic: 2.94 on 5 and 832 DF, p-value: 0.01223
kbl(anova(Elytral), digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria")
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| SeasonDay | 1 | 2.66 | 2.66 | 2.87 | 0.09 |
| Location | 2 | 3.71 | 1.86 | 2.00 | 0.14 |
| SeasonDay:Location | 2 | 7.28 | 3.64 | 3.92 | 0.02 |
| Residuals | 832 | 773.03 | 0.93 | NA | NA |
Use log-transformed elytron length in downstream analysis when using elytron length as an explanatory variable rather than a response variable with respect to variation in body mass due to allometric scaling.
Season Day doesn’t really matter for elytron length (this is expected if elytron length is fixed for the life of adult).
Location doesn’t really matter (caveats: sampling efforts at sites were not equal and non-randomly distributed across the season)
SeasonDay:Location is significant. This is apparent if you look at the violin plots for Season Day grouped by Location. There is a downward trend in BNA, but not in the others. However, the sample size (replicate sampling events at a particular site across the season) is small (BNA = sampled 3 times in the first half of the season, BUF = sampled twice, FDBCC = sampled 4 times, all in second half of season).
Side Note: Vencl and Carlson (1998) (https://www.proquest.com/docview/757065622?accountid=9784) found a bimodal distribution in elytron length at their sites on Long Island across season (July), no effect of year (3 field seasons) so these were all pooled in the distribution displayed in the manuscript. They log-transformed elytron length for their analysis, and the plot displays the untransformed data.
Plots:
#plot with all locations over the season
p_elytral_length_by_season_day <- ggplot(Adult_Final, aes(x=as.numeric(SeasonDay), y=ElytralLength, color = Location)) +
geom_point(size = 1, alpha = 0.7) +
theme_classic() +
ylab("Elytron length (mm)") +
xlab("Season Day") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90), legend.position = "none") +
scale_color_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759"))
p_elytral_length_by_season_day
#BNA
Elytral_BNA<-lm(ElytralLength~SeasonDay, data=Adult_Final[Adult_Final$Location == "USA: Montour Co, Bucknell Natural Area",])
summary(Elytral_BNA)
##
## Call:
## lm(formula = ElytralLength ~ SeasonDay, data = Adult_Final[Adult_Final$Location ==
## "USA: Montour Co, Bucknell Natural Area", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.79806 -0.62992 0.07903 0.66613 2.41613
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.4826 0.1758 53.937 <2e-16 ***
## SeasonDay -0.0243 0.0110 -2.209 0.028 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.954 on 296 degrees of freedom
## Multiple R-squared: 0.01621, Adjusted R-squared: 0.01289
## F-statistic: 4.878 on 1 and 296 DF, p-value: 0.02797
p_elytral_length_by_season_day_BNA <- ggplot(Adult_Final[Adult_Final$Location == "USA: Montour Co, Bucknell Natural Area",], aes(x=as.numeric(SeasonDay), y=ElytralLength)) +
geom_point(size = 1, alpha = 0.7, color = "#4E79A7") +
theme_classic() +
ylab("Elytron length (mm)") +
xlab("Season Day") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90), legend.position = "none") +
geom_smooth(method = "lm", se=FALSE, color = "gray23") +
stat_regline_equation(label.x = 17, label.y = 11.4, aes(label = ..eq.label..)) +
stat_regline_equation(label.x = 17, label.y = 11.0, aes(label = ..rr.label..)) +
ggtitle("BNA")
#p_elytral_length_by_season_day_BNA
Elytral_BUF<-lm(ElytralLength~SeasonDay, data=Adult_Final[Adult_Final$Location == "USA: Union Co, Bucknell Farm",])
summary(Elytral_BUF)
##
## Call:
## lm(formula = ElytralLength ~ SeasonDay, data = Adult_Final[Adult_Final$Location ==
## "USA: Union Co, Bucknell Farm", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.27723 -0.61695 0.02277 0.75613 2.02277
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.41399 0.32250 26.090 <2e-16 ***
## SeasonDay 0.02640 0.01594 1.656 0.102
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.065 on 68 degrees of freedom
## Multiple R-squared: 0.03877, Adjusted R-squared: 0.02463
## F-statistic: 2.743 on 1 and 68 DF, p-value: 0.1023
p_elytral_length_by_season_day_BUF <- ggplot(Adult_Final[Adult_Final$Location == "USA: Union Co, Bucknell Farm",], aes(x=as.numeric(SeasonDay), y=ElytralLength)) +
geom_point(size = 1, alpha = 0.7, color = "#F28E2B") +
theme_classic() +
ylab("Elytron length (mm)") +
xlab("Season Day") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90), legend.position = "none") +
ggtitle("BUF")
#p_log_elytral_length_by_season_day_BUF
Elytral_FDBCC<-lm(ElytralLength~SeasonDay, data=Adult_Final[Adult_Final$Location == "USA: Union Co, Bucknell Ropes Course",])
summary(Elytral_FDBCC)
##
## Call:
## lm(formula = ElytralLength ~ SeasonDay, data = Adult_Final[Adult_Final$Location ==
## "USA: Union Co, Bucknell Ropes Course", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.98747 -0.63048 -0.03074 0.62789 2.91253
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.995642 0.151183 59.502 <2e-16 ***
## SeasonDay -0.002245 0.006267 -0.358 0.72
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9546 on 468 degrees of freedom
## Multiple R-squared: 0.0002742, Adjusted R-squared: -0.001862
## F-statistic: 0.1283 on 1 and 468 DF, p-value: 0.7203
p_elytral_length_by_season_day_FDBCC <- ggplot(Adult_Final[Adult_Final$Location == "USA: Union Co, Bucknell Ropes Course",], aes(x=as.numeric(SeasonDay), y=ElytralLength)) +
geom_point(size = 1, alpha = 0.7, color = "#E15759") +
theme_classic() +
ylab("Elytron length (mm)") +
xlab("Season Day") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90), legend.position = "none") +
ggtitle("FDBCC")
#p_log_elytral_length_by_season_day_FDBCC
p_elytral_length_by_season_day_BNA + p_elytral_length_by_season_day_BUF + p_elytral_length_by_season_day_FDBCC
Conclusions:
#generate distribution of raw data
z <- ggplot(Adult_Final, aes(x=Mass)) +
geom_histogram() +
xlab("Body mass (g)") +
ylab("N fireflies") +
theme(plot.title = element_text(hjust = 0.5))
#generate log distribution plot
y <- ggplot(Adult_Final, aes(x=log(Mass))) +
geom_histogram() +
xlab("log(Body mass) (g)") +
ylab("N fireflies") +
theme(plot.title = element_text(hjust = 0.5))
grid.arrange(z, y, nrow= 1)
#generate tibble
t<- Adult_Final %>%
group_by(Location) %>%
summarise(N=n(), mean_mass=round(mean(Mass),4), SD = round(sd(Mass), 4))
#generate nice looking table
kbl(t, col.names = c("Location", "n", "Mean body mass (g)", "SD"), align = "lccc") %>%
kable_classic(full_width = F, html_font = "Cambria")
| Location | n | Mean body mass (g) | SD |
|---|---|---|---|
| USA: Montour Co, Bucknell Natural Area | 298 | 0.0323 | 0.0092 |
| USA: Union Co, Bucknell Farm | 70 | 0.0326 | 0.0088 |
| USA: Union Co, Bucknell Ropes Course | 470 | 0.0312 | 0.0083 |
#generate sample size labels
site_labels = c(paste("BNA", " (N=", t$N[1], ")", sep=""),
paste0("BUF", " (N=", t$N[2], ")"),
paste0("FDBCC", " (N=", t$N[3], ")"))
#generate text placeholder
i <- ggplot(Adult_Final, aes(x=log(Mass), fill=Location)) +
geom_histogram(alpha= 0.5, position="identity") +
scale_fill_manual(name = "Collection Location", labels = site_labels, values = site_colors)
legend.i <- cowplot::get_legend(i)
#generate histogram with log mass
j <- ggplot(Adult_Final, aes(x=log(Mass), fill=Location)) +
geom_histogram(alpha= 0.5, position="identity") +
xlab("log(Body mass) (g)") +
ylab("N fireflies") +
theme(plot.title = element_text(hjust = 0.5), legend.position = "none") +
scale_fill_manual(name = "Location", labels = site_labels, values = site_colors)
#generate violin plot with log mass
k <- ggplot(Adult_Final, aes(x=Location, y=log(Mass), fill=Location)) +
geom_violin(alpha = 0.9) +
theme_classic() +
ylab("log(Body mass) (g)") +
theme(plot.title = element_text(hjust = 0.5), legend.position = "none", axis.text.x=element_blank(), axis.title.x=element_blank()) +
geom_boxplot(width=0.2, alpha= 0.5, fill = "white") +
scale_fill_manual(values = site_colors)
#generate boxplot
l <- ggplot(Adult_Final, aes(x=Location, y=log(Mass), fill=Location)) +
geom_boxplot(alpha = 0.9) +
geom_jitter(size = 0.1)+
theme_classic() +
ylab("log(Body mass) (g)") +
theme(plot.title = element_text(hjust = 0.5), legend.position = "none", axis.text.x=element_blank(), axis.title.x=element_blank()) +
scale_fill_manual(values = site_colors)
grid.arrange(legend.i,j,k,l, nrow=2)
#generate tibble
t <- Adult_Final %>%
group_by(Location, SeasonDay) %>%
summarise(N=n(), mean_mass=round(mean(Mass),4), SD = round(sd(Mass),4)) %>%
arrange(SeasonDay)
#make table nice in kable
kbl(t, col.names = c("Location", "Season Day", "n", "Mean body mass (g)", "SD"), align = "lcccc") %>%
kable_classic(full_width = F, html_font = "Cambria")
| Location | Season Day | n | Mean body mass (g) | SD |
|---|---|---|---|---|
| USA: Montour Co, Bucknell Natural Area | 9 | 92 | 0.0337 | 0.0092 |
| USA: Union Co, Bucknell Farm | 13 | 47 | 0.0325 | 0.0095 |
| USA: Montour Co, Bucknell Natural Area | 15 | 120 | 0.0313 | 0.0099 |
| USA: Union Co, Bucknell Ropes Course | 17 | 231 | 0.0321 | 0.0087 |
| USA: Montour Co, Bucknell Natural Area | 22 | 86 | 0.0320 | 0.0081 |
| USA: Union Co, Bucknell Ropes Course | 23 | 95 | 0.0312 | 0.0088 |
| USA: Union Co, Bucknell Farm | 30 | 23 | 0.0330 | 0.0074 |
| USA: Union Co, Bucknell Ropes Course | 31 | 90 | 0.0292 | 0.0071 |
| USA: Union Co, Bucknell Ropes Course | 36 | 54 | 0.0306 | 0.0072 |
#generate labels for seasondays that include sample size
t$SeasonDayLabel = paste(t$SeasonDay, " (", t$N, ")", sep ="")
seasonday_label <- c(
'9' = t$SeasonDayLabel[1],
'13' = t$SeasonDayLabel[2],
'15' = t$SeasonDayLabel[3],
'17' = t$SeasonDayLabel[4],
'22' = t$SeasonDayLabel[5],
'23' = t$SeasonDayLabel[6],
'30' = t$SeasonDayLabel[7],
'31' = t$SeasonDayLabel[8],
'36' = t$SeasonDayLabel[9]
)
#generate histogram
m <- ggplot(Adult_Final, aes(x=log(Mass), fill=Location)) +
geom_histogram() +
xlab("log(Body mass) (g)") +
ylab("N fireflies") +
theme(plot.title = element_text(hjust = 0.5)) +
facet_wrap(~SeasonDay, labeller = as_labeller(seasonday_label)) +
scale_fill_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759"))
#generate legend
legend.m <- cowplot::get_legend(m)
#generate histogram
n <- ggplot(Adult_Final, aes(x=log(Mass), fill=Location)) +
geom_histogram() +
xlab("log(Body mass) (g)") +
ylab("N fireflies") +
ggtitle("Season Day (n)") +
theme(plot.title = element_text(hjust = 0.5), legend.position = "none") +
facet_wrap(~SeasonDay, labeller = as_labeller(seasonday_label), ncol = 5) +
scale_fill_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759"))
#generate violin plot by seasonday
o <- ggplot(Adult_Final, aes(x=as.factor(SeasonDay), y=log(Mass), fill = Location)) +
geom_violin(alpha = 0.9) +
theme_classic() +
ylab("log(Body mass) (g)") +
xlab("Season Day") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90), legend.position = "none") +
geom_boxplot(width=0.2, fill = "white") +
scale_fill_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759"))
#generate violin plot by site
p <- ggplot(Adult_Final, aes(x=as.factor(SeasonDay), y=log(Mass), fill = Location)) +
geom_violin(alpha = 0.9) +
theme_classic() +
ylab("log(Body mass) (g)") +
xlab("Season Day") +
theme(plot.title = element_text(hjust = 0.5), legend.position = "none", axis.text.x = element_text(angle = 90)) +
facet_wrap(~Location, labeller = labeller(Location =
c("USA: Montour Co, Bucknell Natural Area" = "BNA",
"USA: Union Co, Bucknell Farm" = "BUF",
"USA: Union Co, Bucknell Ropes Course" = "FDBCC"))) +
geom_boxplot(width=0.2, fill = "white") +
scale_fill_manual(values = c("#4E79A7", "#F28E2B", "#E15759"))
lay1 <- rbind(c(2,2,2,1))
lay2 <- rbind(c(1,1,2,2,2))
g3 <- arrangeGrob(legend.m, n, layout_matrix = lay1)
g4 <- arrangeGrob(o,p, layout_matrix = lay2)
grid.arrange(g3, g4, nrow=2)
Testing the raw data:
#qq plots to assess normality with raw data
q <- ggplot(Adult_Final, aes(sample=Mass, color=Location)) +
stat_qq_band(bandType = "pointwise", fill = "#8DA0CB", alpha = 0.4) +
stat_qq_line(color = "#8DA0CB") +
stat_qq_point() +
facet_wrap(~SeasonDay) +
scale_color_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759")) +
theme(legend.position = "none")
#statistically test
shapiro.table <- NULL
for (day in levels(factor(Adult_Final$SeasonDay))){
#subset table
table.sub <- Adult_Final[Adult_Final$SeasonDay == day,]
test.result <- (shapiro.test(table.sub$Mass))
info <- data.frame(SeasonDay = day, W = round(test.result[1]$statistic, 2), p = round(test.result[2]$p.value, 2))
shapiro.table <- rbind(shapiro.table, info)
}
#generate kable filename
kable_filename = paste("Body_mass_raw_shapiro", ".png", sep="")
#generate nice table with kable
kbl(shapiro.table, row.names = FALSE, col.names = c("Season Day", "W", "P"), align = "ccc") %>%
kable_classic(full_width = F, html_font = "Cambria") %>%
save_kable(file = kable_filename, bs_theme = "cerulean", zoom = 7)
#read in kable table as image
r <- ggdraw() +
draw_image(kable_filename)
#stack the legend on top of the table
g3 <- arrangeGrob(legend.e, r, nrow =2)
#make the column width layout
shapiro_layout <- rbind(c(1,1,1,2))
#plot the graph next to the legend and table
grid.arrange(q, g3, nrow = 1, layout_matrix = shapiro_layout)
Conclusions: Not normal (5/9 days):
However, sample size is large on these days - only Day 30 (not included in list of collections that violate normality) has sample size < 30. Next, test log-transformed data to see if improves.
Note: Log-transformation of body mass is necessary for the assumption of linearity in survival models (see Methods in main text)
Testing log-transformed data:
#qq plots to assess normality
q <- ggplot(Adult_Final, aes(sample=log(Mass), color=Location)) +
stat_qq_band(bandType = "pointwise", fill = "#8DA0CB", alpha = 0.4) +
stat_qq_line(color = "#8DA0CB") +
stat_qq_point() +
facet_wrap(~SeasonDay) +
scale_color_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759")) +
theme(legend.position = "none")
#statistically test
shapiro.table <- NULL
for (day in levels(factor(Adult_Final$SeasonDay))){
#subset table
table.sub <- Adult_Final[Adult_Final$SeasonDay == day,]
test.result <- (shapiro.test(log(table.sub$Mass)))
info <- data.frame(SeasonDay = day, W = round(test.result[1]$statistic, 2), p = round(test.result[2]$p.value, 2))
shapiro.table <- rbind(shapiro.table, info)
}
#generate kable table filename
kable_filename = paste("Body_mass_log_shapiro", ".png", sep="")
#generate nice table with kable
kbl(shapiro.table, row.names = FALSE, col.names = c("Season Day", "W", "P"), align = "ccc") %>%
kable_classic(full_width = F, html_font = "Cambria") %>%
save_kable(file = kable_filename, bs_theme = "cerulean", zoom = 7)
#read in kable table as image
r <- ggdraw() +
draw_image(kable_filename)
#stack the legend on top of the table
g3 <- arrangeGrob(legend.e, r, nrow =2)
#make the column width layout
shapiro_layout <- rbind(c(1,1,1,2))
#plot the graph next to the legend and table
grid.arrange(q, g3, nrow = 1, layout_matrix = shapiro_layout)
Conclusions:
#no Elytral length; + interaction
Mass<-lm(log(Mass)~SeasonDay*Location, data=Adult_Final)
summary(Mass)
##
## Call:
## lm(formula = log(Mass) ~ SeasonDay * Location, data = Adult_Final)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.74613 -0.18636 0.00805 0.18811 0.72058
##
## Coefficients:
## Estimate Std. Error
## (Intercept) -3.4227082 0.0501670
## SeasonDay -0.0034295 0.0031398
## LocationUSA: Union Co, Bucknell Farm -0.0684593 0.0964965
## LocationUSA: Union Co, Bucknell Ropes Course 0.0129793 0.0661484
## SeasonDay:LocationUSA: Union Co, Bucknell Farm 0.0052895 0.0051443
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course -0.0005526 0.0036129
## t value Pr(>|t|)
## (Intercept) -68.226 <2e-16 ***
## SeasonDay -1.092 0.275
## LocationUSA: Union Co, Bucknell Farm -0.709 0.478
## LocationUSA: Union Co, Bucknell Ropes Course 0.196 0.844
## SeasonDay:LocationUSA: Union Co, Bucknell Farm 1.028 0.304
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course -0.153 0.878
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2722 on 832 degrees of freedom
## Multiple R-squared: 0.01091, Adjusted R-squared: 0.00497
## F-statistic: 1.836 on 5 and 832 DF, p-value: 0.1033
#anova(Mass) #interaction term NS, SeasonDay p = 0.008244
#display results as nice table with kable
kbl(anova(Mass), digits = 2) %>%
kable_classic(full_width = F, html_font = "Cambria")
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| SeasonDay | 1 | 0.49 | 0.49 | 6.66 | 0.01 |
| Location | 2 | 0.06 | 0.03 | 0.39 | 0.67 |
| SeasonDay:Location | 2 | 0.13 | 0.06 | 0.87 | 0.42 |
| Residuals | 832 | 61.66 | 0.07 | NA | NA |
Use log-transformed mass in downstream analysis
Mass decreases over the season (potentially reflects using up energy reserves, desiccation).
Plot for multi-panel figure
p_body_mass_by_season_day <- ggplot(Adult_Final, aes(x=as.numeric(SeasonDay), y=log(Mass))) +
geom_point(aes(color = Location), size = 1, alpha = 0.7) +
geom_smooth(method = "lm", se=FALSE, color = "gray23") +
theme_classic() +
ylab("log(Body mass) (g)") +
xlab("Season Day") +
theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90), legend.position = "none") +
scale_color_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759")) + stat_regline_equation(label.x = 23, label.y = -2.8, aes(label = ..eq.label..)) +
stat_regline_equation(label.x = 23, label.y = -2.6, aes(label = ..rr.label..))
p_body_mass_by_season_day
Both<-lm(log(Mass)~log(ElytralLength), data=Adult_Final)
summary(Both)
##
## Call:
## lm(formula = log(Mass) ~ log(ElytralLength), data = Adult_Final)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.86176 -0.13766 0.00558 0.13144 0.66330
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.28523 0.13514 -53.91 <2e-16 ***
## log(ElytralLength) 1.73265 0.06159 28.13 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1957 on 836 degrees of freedom
## Multiple R-squared: 0.4863, Adjusted R-squared: 0.4857
## F-statistic: 791.4 on 1 and 836 DF, p-value: < 2.2e-16
kbl(anova(Both), row.names = TRUE, digits = 130) %>%
kable_classic(full_width = F, html_font = "Cambria")
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| log(ElytralLength) | 1 | 30.31656 | 30.31656157 | 791.4157 | 4.705031e-123 |
| Residuals | 836 | 32.02444 | 0.03830675 | NA | NA |
#using ggpubr package to get equations on graphs
#SOURCE: https://www.roelpeters.be/how-to-add-a-regression-equation-and-r-squared-in-ggplot2/
p_body_mass_by_elytron_length <- ggplot(Adult_Final,aes(x = log(ElytralLength), y = log(Mass))) +
geom_point(data = Adult_Final, aes(color = Location), size = 1, alpha = 0.7) +
geom_smooth(method = "lm", se=FALSE, color = "gray23") +
theme_classic() +
stat_regline_equation(label.y = -2.75, aes(label = ..eq.label..)) +
stat_regline_equation(label.y = -2.87, aes(label = ..rr.label..)) +
scale_color_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759")) +
theme(legend.position = "none") +
ylab("log(Body mass) (g)") +
xlab("log(Elytron length) (mm)")
p_body_mass_by_elytron_length
#including interaction, but restricting to 2nd order interactions due to small sample size
Mass<-lm(log(Mass)~log(ElytralLength)*SeasonDay*Location, data=Adult_Final)
summary(Mass)
##
## Call:
## lm(formula = log(Mass) ~ log(ElytralLength) * SeasonDay * Location,
## data = Adult_Final)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.81332 -0.13162 0.00501 0.12579 0.65748
##
## Coefficients:
## Estimate
## (Intercept) -7.231751
## log(ElytralLength) 1.698810
## SeasonDay -0.032489
## LocationUSA: Union Co, Bucknell Farm 1.997806
## LocationUSA: Union Co, Bucknell Ropes Course -1.037070
## log(ElytralLength):SeasonDay 0.015145
## log(ElytralLength):LocationUSA: Union Co, Bucknell Farm -0.866447
## log(ElytralLength):LocationUSA: Union Co, Bucknell Ropes Course 0.520282
## SeasonDay:LocationUSA: Union Co, Bucknell Farm -0.013209
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course 0.073983
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Farm 0.004831
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course -0.035764
## Std. Error
## (Intercept) 0.848913
## log(ElytralLength) 0.379069
## SeasonDay 0.056310
## LocationUSA: Union Co, Bucknell Farm 1.350131
## LocationUSA: Union Co, Bucknell Ropes Course 1.054078
## log(ElytralLength):SeasonDay 0.025214
## log(ElytralLength):LocationUSA: Union Co, Bucknell Farm 0.612366
## log(ElytralLength):LocationUSA: Union Co, Bucknell Ropes Course 0.474623
## SeasonDay:LocationUSA: Union Co, Bucknell Farm 0.078929
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course 0.061974
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Farm 0.035582
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course 0.027855
## t value
## (Intercept) -8.519
## log(ElytralLength) 4.482
## SeasonDay -0.577
## LocationUSA: Union Co, Bucknell Farm 1.480
## LocationUSA: Union Co, Bucknell Ropes Course -0.984
## log(ElytralLength):SeasonDay 0.601
## log(ElytralLength):LocationUSA: Union Co, Bucknell Farm -1.415
## log(ElytralLength):LocationUSA: Union Co, Bucknell Ropes Course 1.096
## SeasonDay:LocationUSA: Union Co, Bucknell Farm -0.167
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course 1.194
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Farm 0.136
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course -1.284
## Pr(>|t|)
## (Intercept) < 2e-16
## log(ElytralLength) 8.46e-06
## SeasonDay 0.564
## LocationUSA: Union Co, Bucknell Farm 0.139
## LocationUSA: Union Co, Bucknell Ropes Course 0.325
## log(ElytralLength):SeasonDay 0.548
## log(ElytralLength):LocationUSA: Union Co, Bucknell Farm 0.157
## log(ElytralLength):LocationUSA: Union Co, Bucknell Ropes Course 0.273
## SeasonDay:LocationUSA: Union Co, Bucknell Farm 0.867
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course 0.233
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Farm 0.892
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course 0.200
##
## (Intercept) ***
## log(ElytralLength) ***
## SeasonDay
## LocationUSA: Union Co, Bucknell Farm
## LocationUSA: Union Co, Bucknell Ropes Course
## log(ElytralLength):SeasonDay
## log(ElytralLength):LocationUSA: Union Co, Bucknell Farm
## log(ElytralLength):LocationUSA: Union Co, Bucknell Ropes Course
## SeasonDay:LocationUSA: Union Co, Bucknell Farm
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Farm
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1933 on 826 degrees of freedom
## Multiple R-squared: 0.5048, Adjusted R-squared: 0.4982
## F-statistic: 76.53 on 11 and 826 DF, p-value: < 2.2e-16
kbl(anova(Mass), row.names = TRUE, digits = 130) %>%
kable_classic(full_width = F, html_font = "Cambria")
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| log(ElytralLength) | 1 | 30.31656157 | 30.31656157 | 811.0826030 | 7.886660e-125 |
| SeasonDay | 1 | 0.16518932 | 0.16518932 | 4.4194385 | 3.583390e-02 |
| Location | 2 | 0.29247175 | 0.14623587 | 3.9123623 | 2.036480e-02 |
| log(ElytralLength):SeasonDay | 1 | 0.02514312 | 0.02514312 | 0.6726735 | 4.123577e-01 |
| log(ElytralLength):Location | 2 | 0.39495934 | 0.19747967 | 5.2833276 | 5.248466e-03 |
| SeasonDay:Location | 2 | 0.15246708 | 0.07623354 | 2.0395354 | 1.307438e-01 |
| log(ElytralLength):SeasonDay:Location | 2 | 0.12006939 | 0.06003470 | 1.6061550 | 2.012837e-01 |
| Residuals | 826 | 30.87414249 | 0.03737790 | NA | NA |
log(Elytron length) sig. correlated with body mass, but can’t explain everything (v small p, r^2 = 0.49)
Location and SeasonDay significant when elytra length accounted for
Question: take residuals as “adult reserves” measurement?
#using patchwork
#get the plots together
multiplot <- p_elytral_length_by_season_day_BNA + p_elytral_length_by_season_day_BUF + p_elytral_length_by_season_day_FDBCC + p_body_mass_by_season_day + p_body_mass_by_elytron_length
#add annotations and layouts
w <- multiplot + plot_annotation(tag_levels = 'A') + plot_layout(widths = unit(c(3),c("in")), heights = unit(c(2), c("in")), guides = "collect") & theme(legend.position = 'bottom')
w
#ggsave("FigS1.EL_Mass_analysis.png", w, width = 12, height = 8, device = "png")
sessionInfo()
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] patchwork_1.1.2 ggpubr_0.4.0 cowplot_1.1.1 sjPlot_2.8.11
## [5] kableExtra_1.3.4 gridExtra_2.3 fabricatr_1.0.0 ggthemes_4.2.4
## [9] qqplotr_0.0.5 Hmisc_4.7-1 Formula_1.2-4 survival_3.4-0
## [13] lattice_0.20-45 dplyr_1.0.10 tidyr_1.2.1 ggplot2_3.3.6
##
## loaded via a namespace (and not attached):
## [1] minqa_1.2.4 colorspace_2.0-3 ggsignif_0.6.3
## [4] deldir_1.0-6 ellipsis_0.3.2 sjlabelled_1.2.0
## [7] estimability_1.4.1 htmlTable_2.4.1 parameters_0.18.2
## [10] base64enc_0.1-3 rstudioapi_0.14 farver_2.1.1
## [13] fansi_1.0.3 mvtnorm_1.1-3 xml2_1.3.3
## [16] splines_4.2.1 cachem_1.0.6 robustbase_0.95-0
## [19] knitr_1.40 sjmisc_2.8.9 polynom_1.4-1
## [22] jsonlite_1.8.0 nloptr_2.0.3 ggeffects_1.1.3
## [25] broom_1.0.1 cluster_2.1.4 png_0.1-7
## [28] effectsize_0.7.0.5 compiler_4.2.1 httr_1.4.4
## [31] sjstats_0.18.1 emmeans_1.8.1-1 backports_1.4.1
## [34] Matrix_1.5-1 fastmap_1.1.0 cli_3.4.0
## [37] htmltools_0.5.3 tools_4.2.1 gtable_0.3.1
## [40] glue_1.6.2 Rcpp_1.0.9 carData_3.0-5
## [43] jquerylib_0.1.4 vctrs_0.4.1 svglite_2.1.0
## [46] nlme_3.1-159 insight_0.18.4 xfun_0.33
## [49] stringr_1.4.1 ps_1.7.1 lme4_1.1-30
## [52] rvest_1.0.3 lifecycle_1.0.2 rstatix_0.7.0
## [55] DEoptimR_1.0-11 MASS_7.3-58.1 scales_1.2.1
## [58] RColorBrewer_1.1-3 yaml_2.3.5 sass_0.4.2
## [61] rpart_4.1.16 latticeExtra_0.6-30 stringi_1.7.8
## [64] highr_0.9 bayestestR_0.13.0 checkmate_2.1.0
## [67] boot_1.3-28 rlang_1.0.5 pkgconfig_2.0.3
## [70] systemfonts_1.0.4 evaluate_0.16 purrr_0.3.4
## [73] htmlwidgets_1.5.4 labeling_0.4.2 tidyselect_1.1.2
## [76] processx_3.7.0 magrittr_2.0.3 R6_2.5.1
## [79] magick_2.7.3 generics_0.1.3 pillar_1.8.1
## [82] foreign_0.8-82 withr_2.5.0 mgcv_1.8-40
## [85] datawizard_0.6.0 abind_1.4-5 nnet_7.3-17
## [88] tibble_3.1.8 performance_0.9.2 modelr_0.1.9
## [91] car_3.1-0 interp_1.1-3 utf8_1.2.2
## [94] rmarkdown_2.16 jpeg_0.1-9 data.table_1.14.2
## [97] callr_3.7.2 digest_0.6.29 webshot_0.5.3
## [100] xtable_1.8-4 munsell_0.5.0 viridisLite_0.4.1
## [103] bslib_0.4.0